Mon. Not. R. Astron. Soc. , 1-?? (2002) Printed 2 February 2008 (MN WT&L style file v2.2) 



The local sub-mm luminosity functions and 
predictions from Spitzer to Herschel 

Stephen Serjeant 1 and Diana Harrison 2 ' 3 

1 Centre for Astrophysics and Planetary Science, School of Physical Sciences, University of Kent, Canterbury, Kent, CT2 7NZ, UK 

2 Astrophysics Group, Blackett Laboratory, Imperial College, Prince Consort Road, London SW7 2BW, UK 

3 Dept. of Physics, Cavendish Laboratory, Madingley Road, Cambridge CB3 OHE, UK 



Accepted; Received; in original form 2003 March 4 



ABSTRACT 

We present new determinations of the local sub-mm luminosity functions, 
solving the "sub-mm Olbers' Paradox." We also present predictions of 
source counts and luminosity functions in current and future far-infrared 
to sub-mm surveys. Using the sub-mm colour temperature relations from 
the SCUBA Local Universe Galaxy Survey, and the discovery of excess 
450/im excess emission in these galaxies, we interpolate and extrapolate 
the IRAS detections to make predictions of the SEDs of all 15411 PSC-z 
galaxies from 50 — 1300/mi. Despite the long extrapolations we find excel- 
lent agreement with (a) the 90/xm luminosity function of Serjeant et al. 
(2001), (b) the 850/xm luminosity function of Dunne et al. (2000), (c) the 
mm- wave photometry of Andreani & Franceschini (1996); (d) the asymp- 
totic differential and integral source count predictions at 50 — 1300/im 
by Rowan-Robinson (2001). We find the local 850/im sub-mm luminosity 
density converges to 7.3 ±0.2 x 10 19 W Hz -1 Mpc -3 . Remarkably, the 
local spectral luminosity density and the extragalactic background light 
together strongly constrain the cosmic star formation history for a wide 
class of evolutionary assumptions. We find that the extragalactic back- 
ground light, the 850/im 8mJy source counts, and the r2* constraints all 
independently point to a decline in the comoving star formation rate at 
z > 1. In order to reconcile this with direct determinations, we suggest 
either there is a top-heavy initial mass function at high redshifts, and/or 
there is stronger evolution in the more luminous far-infrared galaxies than 
seen in the population as a whole. 

Key words: cosmology: observations - galaxies: formation - infrared: 
galaxies - galaxies: evolution - galaxies: starburst - galaxies: statistics 



1 INTRODUCTION 

The high expectations for sub-mm/mm-wave surveys 
have been fully and variously realised. The discovery 
of a population of z > 1 sub-mm galaxies resolving up 
to 50% of the extragalactic background light, and the 
surprising lack of overlap with the hard X-ray Chan- 
dra populations, have been among the key results in 
cosmology of the past five years. The SCUBA survey 
of the Hubble Deep Field (HDF) demonstrated the 
feasibility of blank-field surveys (Hughes, Serjeant, 
Dunlop, Rowan-Robinson, et al. 1998, Serjeant et 



al. 2003), and several complementary sub-mm/mm- 
wave surveys have since been conducted to a vari- 
ety of depths and areal coverages (e.g. Scott et al. 
2002; Fox et al. 2002; Eales et al. 2000; Barger et al. 
1999). These surveys will be augmented with much 
larger surveys with (e.g.) SCUBA-2 (e.g. Holland et 
al. 2002), BOLOCAM on the CSO (e.g. Glenn et 
al. 1998), the Balloon-borne Large Aperture Sub-mm 
Telescope (e.g. Tucker et al. 2004), the Large Millime- 
ter Telescope (e.g. Kaercher et al. 2000), the Atacama 
Large Millimeter Array (ALMA), Planck High Fre- 
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quency Instrument (HFI, e.g. Lamarre et al. 2001), 
and Herschel SPIRE (e.g. Griffin et al. 2001). 

It is not widely appreciated that the most rea- 
sonable prior expectation for the SCUBA-HDF sur- 
vey, based on the observed optical HDF population, 
was zero sources at 850/im: for example, Thomp- 
son et al. (2001) model the near-IR/optical/near-UV 
SEDs of the NICMOS-HDF galaxies, and derive sub- 
mjy/ micro- Jy fluxes at 850/im. The mjy-level sub- 
mm galaxies found in blank-field surveys on the whole 
can be reasonably regarded as a new galaxy popula- 
tion. 

One attractive interpretation of this new pop- 
ulation is that it is comprised at least in part of 
proto-spheroids (e.g. Dunlop 2001): the star forma- 
tion rates inferred in sub-mm survey galaxies (iS 1000 
M Q /yr) are sufficient to assemble a giant elliptical 
galaxy in < lGyr. In support of this interpreta- 
tion, the K-band morphologies of at least some re- 
semble high-redshift radiogalaxies (Lutz et al. 2002), 
and provisional measurements of the clustering are 
consistent with a high bias parameter population 
at high redshift (Scott et al. 2002). Nevertheless, 
the sub-mm galaxies are not unambiguously identi- 
fied with proto-spheroids; Rowan- Robinson (2001), 
Efstathiou & Rowan-Robinson (2002) and King & 
Rowan-Robinson (2003) instead postulate that the 
blank-field sub-mm survey galaxies are dominated in 
the sub-mm by cirrus heated by their interstellar ra- 
diation fields, rather than ultraluminous star forma- 
tion. In support of this interpretation, Efstathiou & 
Rowan-Robinson (2002) model the SEDs of sub-mm 
galaxies and find them consistent with their cirrus ra- 
diative transfer models. The existence of cool colour 
temperature components in local galaxies, physically 
distinct from the warmer dust dominating the 60/im 
emission, have been pointed out by many authors 
(e.g. Dunne et al. 2000 and refs. therein, Stickel et al. 
2000, Siebenmorgen et al. 1999, Domingue et al. 1999, 
Trewhella et al. 2000; Haas et al. 2000.) If this cool 
component dominates the observed sub-mm fluxes of 
high-z blank field sub-mm/mm-wave survey sources, 
then the emission should be extended over ~ 1" 
rather than confined to a compact ~ 0.1" ultralumi- 
nous starburst (Efstathiou & Rowan- Robinson 2002, 
Lee & Rowan- Robinson 2003, Kaviani et al. 2002). 
This distinction will be possible with the planned 
ALMA. 

In the meantime, the evolution of the sub-mm 
galaxy population can be strongly constrained by the 
integrated extragalactic background light, the local 
multiwavelength luminosity functions, and the source 
counts. The local 850/im luminosity function was de- 
rived in the SCUBA Local Universe Galaxy Survey 
(SLUGS, Dunne et al. 2000) from their SCUBA pho- 
tometry of the IRAS Bright Galaxy Survey. A curious 
aspect of their luminosity function was that the faint 
end slope was not sufficiently shallow for the local 
luminosity density to converge, which the authors re- 



ferred to as the "sub-mm Olbers' paradox". This is 
a pity from the point of view of modelling the high 
redshift population, since the integrated extragalac- 
tic background light is a key constraint. In order to 
find the expected flattening of the luminosity func- 
tion slope at lower luminosities, the SLUGS survey 
is currently being extended with SCUBA photometry 
of optically-selected galaxies. Meanwhile, several au- 
thors have attempted to use the colour temperature 
- luminosity relation found in SLUGS to transform 
the 60/im luminosity function to other wavelengths, 
and hence constrain the high-redshift evolution (e.g. 
Lagache, Dole & Puget 2003, Lee & Rowan- Robinson 
2003, Chapman et al. 2003). The discovery of an ad- 
ditional excess component at 450/im (Dunne & Eales 
2001) relative to their initial colour temperature - lu- 
minosity relation is only rarely included in such mod- 
els (e.g. Lagache, Dole & Puget 2003), and there are 
problems in accurately representing the population 
mix of galaxies. 

In this paper we take an alternative approach 
to determining the multiwavelength local luminosity 
functions. We model the spectral energy distributions 
(SEDs) of all 15411 PSC-z galaxies (Saunders et al. 
2000), constrained by all available far- infrared and 
sub-mm colour-colour relations from SLUGS and else- 
where. This guarantees the correct local population 
mix at every wavelength and minimises the assump- 
tions about the trends of SED shape with luminos- 
ity (see e.g. Rowan- Robinson et al. 2001, in which 
the differences in assumed population mix result in 
rather different conclusions to this work). Our phe- 
nomenological approach is sufficient to determine (for 
example) the 850/im luminosities to within a factor 
of 2 in individual galaxies. While this is not accu- 
rate enough to serve as calibrations for sub-mm mis- 
sions, except statistically over the whole population, 
it is more than sufficiently accurate for determining 
local luminosity functions, which we present here from 
far-infrared to millimetric wavelengths. We use these 
local luminosity functions to model the high- 2 evo- 
lution, constrained by the extragalactic background 
light spectrum, the total cosmological matter density 
in stars O*, and the sub-mm source counts. 

Our SED models are accurate enough to identify 
potential calibrators for further photometric study. 
The SEDs can also be used to determine the bright 
source contributions to current and future far-infrared 
and sub-mm/mm-wave surveys. Predictions for bright 
source catalogues are also essential tools for the plan- 
ning and execution of future wide area far-infrared 
missions. For example, the catalogues are useful for 
both tertiary flux calibration and astrometric calibra- 
tion in all-sky surveys (e.g. Planck HFI, Astro-F, e.g. 
Pearson et al. 2001, Shibai 2002), and also for mission 
strategy (such as deep survey fields avoiding bright 
sources) . 

Section 2 describes our modelling of the spectral 
energy distributions, including observational tests of 
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our methodology. Section 3 presents the results, such 
as multiwavelength local luminosity functions, and 
bright source counts for Planck HFI, Herschel SPIRE, 
and Astro-F The implications of our results are dis- 
cussed in section 4. 

Throughout this paper we adopt the "concor- 
dance" cosmology of Om = 0.3, Qa = 0.7 and 
H = 65 km s" 1 Mpc -1 . 



2 METHOD 



2.1 850/im predictions 

The PSC-z survey is a redshift survey of IRAS galax- 
ies, to a depth of 0.6Jy at 60/tm. PSC-z covers 84% 
of the sky, and has a median redshift of 8400 km s _1 . 
Further details can be found in Saunders et al. (2000) . 

Our first stage in modelling the far-infrared to 
mm-wave SEDs of IRAS PSC-Z galaxies is to deter- 
mine their 850/im fluxes. To do this we use the re- 
sults from the 850/im photometry programme of the 
IRAS Bright Galaxy Sample by Dunne et al. (2000). 
These authors found a strong correlation between far- 
infrared luminosity and sub-mm:far-infrared colour. 
The far-infrared luminosity correlates strongly with 
the 60 — 100/im colour, so an alternative approach to 
modelling the luminosity-colour relation is to fit to 
the 100 — 60/im vs. 850 — 60/im colour-colour plane. 
This is plotted in figure 1, together with the predic- 
tions from grey-body emission: 



S v = v B B v = v 1 



,2/iz/ 3 



exp(hv/{kT)) - 1 



(1) 



The lines in figure 1 correspond to (3 = 
1.1, 1.2, 1.3, 1.4, 1.5. Dunne et al. (2000) also modelled 
the colour temperatures and /3 values for all their 
sample using limited 450/im data, and found a me- 
dian f3 = 1.3 which is in excellent agreement with the 
best fit to figure 1. Physically, these galaxies are very 
unlikely to be single temperature or even dual tem- 
perature systems (e.g. Efstathiou, Rowan-Robinson 
& Siebenmorgen 2000) but all that matters for our 
current purpose is that we have a phenomenological 
model which can be used to interpolate to other fre- 
quencies. 

We use the co-added or extended ADDSCAN 
IRAS fluxes (Saunders et al. 2000). For the 29% of 
sources with only 100/im limits, we use the non-linear 
60/im-100/im luminosity-luminosity correlation (i.e., 
the well-known trend of IRAS colour temperature 
correlating with luminosity) to estimate the 100/im 
fluxes. Redshift information is available for almost 
all of PSC-Z, but because redshift is degenerate with 
colour temperature this information is not necessary 
for sub-mm flux prediction where 100/im detections 
are available. 



Figure 1. Sub-mm:far-infrarcd colour-colour plane for the 
Dunne et al. (2000) sample. From top to bottom, lines 
show the predictions for grey body indices of (3 = 1.1, 1.2, 
1.3, 1.4, 1.5. 



2.2 450/im predictions and 

interpolations /extrapolations 

While the phenomenological model above is a good 
predictor of the 850/im fluxes of IRAS galaxies, there 
is an additional component at shorter wavelengths 
discovered by Dunne & Eales (2001). 

The combined 850/im, 450/im, 100/im and 60/im 
data of the IRAS Bright Galaxy Survey led Dunne & 
Eales (2001) to argue for two (3 = 2 grey-body compo- 
nents, because the excess flux in the short-wavelength 
sub-mm data could not be fit by any single grey- 
body template. As the two component model is a 
four-parameter fit (two temperatures and two nor- 
malisations) with the emissivity fixed, four-band pho- 
tometry is just sufficient in principle to define the fit 
uniquely. 

Together with the IRAS measurements, our pre- 
scription for 850/xm fluxes gives us a total of three 
bands on which to base SED interpolations. We can 
obtain a further phenomenological prediction for the 
450/im fluxes on the linear 60 : 850/im vs. 60 : 450/im 
colour-colour correlation, for which the IRAS Bright 
Galaxy Survey yields 



log 10 (S 6 o/S 4 5o) = log 10 (S , 6o/S , 85o) - 0.924 



(2) 



These prescriptions should be sufficient to define the 
sub-mm fluxes to within a factor of two. 

We then numerically solve for the parameters of 
the two temperature components, using the four-band 
measured and predicted photometry. Note that this 
two component fit uses a grey body index of (5 = 2, 
which is not consistent with the f3 = 1.3 model used 
to predict the 850/im fluxes. This does not affect the 
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self-consistency or utility of the two-component fit, 
because the only relevant information is that a numer- 
ical scheme of whatever nature is available to predict 
the 850/wn fluxes of the IRAS PSC-Z galaxies. Also, it 
is worth stressing again that more physically reason- 
able models incorporate a continuum of temperature 
components; all that matters for our current purpose 
is to obtain a phenomenological parametric model on 
which to base predictions of the far-IR to mm-wave 
SEDs. In the following section we will test the con- 
sistency of the phenomenological predictions using a 
variety of observational constraints. 

The dual-temperature model incorporates the 
available photometric constraints on the far-infrared 
to mm-wave SEDs of IR-luminous galaxies, and pro- 
vides the best available interpolations and extrapola- 
tions from 50 — 1300/xm in an all-sky catalogue. The 
PSC-Z predictions for Astro-F, Herschel, Planck HFI 
and ISO are available as an IDL database from the 
authors. Also available is interpolation software to 
provide further predictions at further arbitrary wave- 
lengths. 



3 RESULTS 

3.1 Observational tests 

In figure 2 we plot the observed 1.25mm photometry 
of Andreani & Franceschini (1996), compared with 
our predictions. These observations are subject to a 
systematic uncertainty of up to a factor of 2 due to 
aperture corrections. Not withstanding a large sys- 
tematic shift in flux calibration, our predictions are 
in good agreement with these observations. 

A stronger tests of our methodology comes from 
the Serendipity Survey of the Infrared Space Obser- 
vatory (ISOSS), taken at 175^m by the ISOPHOT 
instrument (e.g. Stickel et al. 2000). The survey cov- 
ered ~ 10% of the sky, and although the completeness 
is not yet well understood the flux calibration is well- 
determined enough to provide a useful test. Figure 3 
shows the predicted 175/nm fluxes against those ob- 
served. The agreement is excellent. This also demon- 
strates that the cool dust excess reported by Stickel 
et al. (2000) is identical to that reported by Dunne & 
Ealcs (2001), at least insofar as fitting of SEDs with 
colour temperature components is concerned. (This 
does not necessarily imply the 175/im-emitting dust 
is cospatial to the 450/im-emitting dust, e.g. section 
3.3.) 

There are few other examples of objects with 
sufficient multi-wavelength coverage to provide tests 
of our methodology (a problem also encountered by 
other authors: e.g. Hughes et al. 2002). In figure 4 
we plot the ratio of our predictions with the obser- 
vations for the SLUGS survey, as a function of wave- 
length. The dispersion is a measure of the error in 
our predictions for individual objects, while the mean 
values demonstrate that the systematic errors in our 



Figure 2. Comparison of predicted and observed 1.25mm 
fluxes for the sample of Andreani & Franceschini 1996, 
using their a mm = a /3 aperture correction. Their alter- 
native a m m = ct aperture corrections result in an upward 
revision of the observed fluxes by a factor of ~ xl.5 — 2. 
The straight line indicates the locus of exact agreement 
between observations and predictions. 

methodology are typically only at the few percent 
level in the wavelength range covered. Our system- 
atic agreement is slightly better than that of Dale et 
al. (2001, 2002) who use a single-component SED but 
with a varying grey body index. 

One of the individual galaxies in the SLUGS sur- 
vey, also covered by the ISOSS survey, is shown in 
figure 5. The agreement is again very good. Note that 
this is not a fit to the multi-wavelength data; instead, 
the model is the mean SED for a galaxy with the 
60 : 100/im colour of this galaxy. 

3.2 Multi-wavelength local luminosity 
functions 

Our sample is 60/im-selected, but that does not pre- 
clude us from constructing a luminosity function at 
another wavelength provided there are no large popu- 
lations that would be missing from our sample at any 
redshift, compared to a flux limited sample at that 
other wavelength. Serjeant et al. (2001) describe the 
methodology for deriving luminosity functions in the 
face of complicated selection functions, and the condi- 
tions for the validity of such derivations. In summary, 
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Figure 3. Comparison of predicted and observed 175/im 
fluxes for the ISOPHOT Serendipity Survey of Stickel et 
al. (2000), using the best available current flux calibration 
(Stickel. priv. comm.). 

the number density in a given luminosity bin is given 
by 

$ = ECKnax^r 1 (3) 

and the error on this quantity is given by 

A$ = VS(\/ ma x,0~ 2 (4) 

Here, Knax.i = V(z max ,i) where V is the comoving 
volume, and z ma x,i is the redshift at which object i 
would be observed at the 60pm PSC-z flux limit, if 
all other intrinsic properties of that object were kept 
the same. Our methodology differs slightly from that 
of Takeuchi et al. (2003), in that those authors used 
different K-corrections and a different sub-samples 
(warm and cool) of the PSC-z survey, as well as mak- 
ing a correction for density fluctuations. These differ- 
ences in methodology lead to small changes in (for 
example) the faint end slope, but are too small to 
affect the conclusions of our paper. 

Although we have found that the 850pm predic- 
tions for individual galaxies are only good to within 
a factor of ~ 2 (figure 4), this is nevertheless accu- 
rate enough for the determination of the local lumi- 
nosity function. Note that these luminosities are esti- 
mated directly from the colour-colour relations, and 
are not interpolated using our two-component dust 
models (except for the K-correction terms). In fig- 



Figure 4. Comparison of the observed SEDs of SLUGS 
galaxies with the predictions of this paper. The X sym- 
bols indicate the SLUGS measurements themselves, while 
the + symbols represent the supplementary data compiled 
from the literature by Dunne & Ealcs (2001). 

ure 6, we plot the l/Umax 850pm luminosity function 
for PSC-Z assuming (1 + z) 3 pure luminosity evolu- 
tion (an assumption which makes a small correction 
to only the brightest data points), and compare it to 
the observed luminosity function from Dunne et al. 
(2000). The agreement is excellent, and as with fur- 
ther multi-wavelength luminosity functions below is 
sufficient to determine the convergence of the local 
luminosity density. This is shown in the case of the 
850pm luminosity density in figure 7, where the PSC- 
z constraints are compared to the existing SLUGS 
data. The SLUGS data reaches the peak of the lu- 
minosity density, but would need to extend around 
~ xlO fainter to determine the convergence in the 
background. Using PSC-z, we find that local luminos- 
ity density at 850pm converges at 7.3 ±0.2 x 10 19 h 6 s 
W Hz -1 Mpc -3 , solving the sub-mm Olbers' Paradox. 
The constraints at other wavelengths will be discussed 
in further detail in section 3.3. 

Unlike the 850pm luminosity function, the 90pm 
luminosities do not have direct colour-colour relations 
to use for their prediction. We therefore resort to 
the two-component SEDs to derive interpolated 90pm 
fluxes. In figure 8 we plot the predicted 90pm lumi- 
nosity function for PSC-Z, and compare it to the ob- 
served luminosity function from the European Large 
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Figure 12. Continuation of figure 11. 
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Table 1. Best-fit parameters for the luminosity functions shown in figure 9. The parameters arc defined in equations 5 
and 6. 



Figure 5. Comparison of the observed SED of NGC520 
with the predictions. The cool and warm components are 
shown separately, and the sum (also shown) is in good 
agreement with the observations. References for the pho- 
tometry are in Dunne & Eales (2001). 

Area ISO Survey (Serjeant et al. 2001). Pure lumi- 
nosity evolution of (1 + z) s was assumed. 

The 175pm luminosity function cannot be deter- 
mined directly from the ISOSS survey, because the 
completeness of this survey is not well-understood. 
In figure 9 we plot the predicted 175/wn luminosity 
function from our PSCZ models, and further exam- 
ples of the interpolated luminosity functions possible 
from our PSC-z SED ensemble. The 175/im luminos- 
ity function is currently the best constraint available 
at around this wavelength before Spitzer and Astro-F, 



Figure 6. Projected local 850/im luminosity function from 
PSC-Z (error bars) assuming pure luminosity evolution of 
(1 + z)' A , compared with the directly measured luminosity 
function from Dunne et al. 2000 (plotted as open circles). 
Also plotted are the best fit Schcchtcr and double power 
law functions (dashed and dotted lines respectively). 



because the number of local 175/im-selected galaxies 
from the European Large Area ISO Survey (Oliver 
et al. 2000) and FIRBACK (Dole et al. 2001) is too 
small for the construction of a luminosity function. 

Table 1 shows the min-x 2 fit parameters for the 
local luminosity functions discussed in this section, 
using the following functional forms: 

$ = $»,s ln(W)(L/L«,s) 1 - a e- L/L ^ s (5) 

$ = $»,P/((L/L*,pf + (L/£,,p) 7 ) (6) 
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Figure 7. Projected local 850/im luminosity density func- 
tion from PSC-Z (error bars), compared with the directly 
measured luminosity density from Dunne et al. 2000 (plot- 
ted as open circles). This is derived directly from figure 6 
after removal of the sr _1 factor in the luminosities. Note 
that the addition of the PSC-z data adds the first evidence 
of convergence in the local luminosity density at 850^tm. 



Figure 8. Projected local 90/^m luminosity function from 
PSC-Z (error bars), compared with the directly measured 
luminosity function from ELAIS (shaded area) from Ser- 
jeant ct al. 2001. In both cases pure luminosity evolution 
of (1 + z) 3 was assumed. Also plotted are the best fit 
Schechter and double power law functions (dashed and 
dotted lines respectively). 



The Schechter functions give acceptable \ 2 values 
when restricted to data less than a few times L* , but 
in no cases provide an acceptable fit at the brightest 
end. The double power law fit, on the other hand, fits 
all the available data. 

3.3 The local extragalactic luminosity 
density 

The extraordinarily tight constraints on the far- 
infrared to sub-mm luminosity functions in section 
3.2 allow us to make tight constraints on the local ex- 
tragalactic spectral luminosity density. This is plot- 
ted in figure 10. The vertical lines indicate where di- 
rect constraints on the PSC-z luminosity densities are 
available via IRAS, SCUBA or ISO. These are spaced 
roughly in factors of 2 in wavelength. Between these 
wavelengths, the luminosity densities are based on 
interpolated flux densities using our two-component 
SED models. The thickness of the line is the ±1<t error 
on the luminosity density. 

Remarkably, the local spectral luminosity density 
itself can be described to within better than a percent 
by two grey bodies, despite being comprised of the 



sum of > 30000 individual grey bodies: 

3.87 x 10 31 
+ A 4 -50( ex p(318.0/A) - 1) { ' 

where A is in pirn, and / is in W/Hz/Mpc 3 . This is an 
instructive lesson: one should not use the two compo- 
nents of our model (or that of Dunne & Eales 2001) 
to argue for two discrete, physically distinct phases. 
The two components are a simple phenomenological 
model, and do not preclude a continuum of dust tem- 
peratures being present. It is also easy to demonstrate 
that the results of radiative transfer models of star 
forming regions (e.g. Efstathiou et al. 2000) can of- 
ten be fit to within a few percent by a small number 
of grey bodies, even though a continuum of tempera- 
tures is intrinsic to the models. 

Finally, one can use equation 7 to calculate the 
local bolometric luminosity density from thermal dust 
emission. In the interval 30/im-3mm, this is found to 
be 2.6 x 10 34 W Mpc~ 3 , i.e. 6.4 x 10 7 L Q Mpc" 3 , 
in very good agreement with the determination by 
Dwek et al. 1998. (Differences in the SED ensemble 
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45.2 
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4.92 



Table 2. Approximate asymptotic Eucludcan source 
count slopes for the predicted source counts plotted in fig- 
ure 11. These asymptotic scounts follow dN/dS = kS~ 2,5 , 
with the coefficient k listed in the table as a function of 
wavelength. 



Figure 10. Local spectral luminosity density. The thick- 
ness of the curve represents the ±1<t error. Vertical lines in- 
dicate the wavelengths where direct estimates are available 
from colour-colour relations using PSC-Z, or wavelengths 
where direct observational tests verify the methodology. 
The values of the local luminosity density at intervening 
wavelengths use the two component dust models of PSC-z. 



do however generate differences in our local spectral 
luminosity density, compared to that of Dwek et al. 
1998.) The SED models however are not constrained 
over all this range; the local luminosity density over 
the range 60- 1300^m is 5.3 x 1O 7 L Mpc~ 3 . Our lo- 
cal bolometric luminosity density differs from that of 
Soifer and Neugebauer (1991) owing to the different 
bolometric corrections for each galaxy, but our inte- 
grated luminosity density over the range 50 — 100/im 
(i.e. the overlap range in measurements and model 
validity) is in good agreement: 4.9 x 1O 6 L0 Mpc~ 3 , 
compared to the result of integrating the Soifer & 
Neugebauer data over this range to obtain 6.3 x 10 6 Lq 
Mpc~ 3 . The integrated background in the optical is 
much larger: 1.26 x W 8 Lq Mpc -3 (Glazebrook et al. 
2003). The fact that the extragalactic background has 
a larger luminous energy fraction in the far-infrared 
implies stronger relative evolution in the far-infrared 
luminous population (e.g. Puget et al. 1996, Hauser et 
al. 1998, Schlegel et al. 1998, Rowan-Robinson 2001, 
Smail et al. 2002). 



3.4 Bright source counts from Astro-F to 
Herschel 

Figures 11 and 12 show the bright source counts 
at wavelengths covered by Spitzer, Astro-F, BLAST 
(Tucker et al. 2004), SCUBA, MAMBO, Herschel, 
and Planck. The bright slopes are Euclidean as ex- 
pected (table 2), but at faint fluxes the counts fall 
below the Euclidean extrapolation. This is an incom- 
pleteness due to the 60/^m flux limit of PSC-z: of the 
local galaxies at the faint end of the sub-mm counts, 
only the warmer will appear in a 60/im-limited sur- 
vey. It is possible to correct for this using accessible- 
volume weightings. This is essentially equivalent to 
constructing luminosity functions, discussed above, 
but the resulting statistic (confirmation of Euclidean 
slope) would have less utility and interest than the 
luminosity functions themselves. Furthermore, PSC-z 
is not deep enough to detect the sub-mm blank-field 
survey population, both because of the less favourable 
K-corrections and because of the low median redshift 
of PSC-z itself. Therefore, deviations from the Eu- 
clidean slope caused by high-z evolution will not be 
well-determined by PSC-z. 

How similar are our PSC-Z "simulated" cata- 
logues at each wavelength to flux-limited surveys at 
those wavelengths? There are two types of object 
which could be undetected by IRAS: local galaxies 
too faint in the IRAS bands due to their intrinsic lu- 
minosities and colour temperatures, and high-redshift 
galaxies beyond the upper redshift limit of PSC-Z. An 
indication of the fluxes at which local galaxies start 
to be missed is the point at which the counts depart 
from the Eucludean slope, in figures 11 and 12. Note 
that this is not an effect which would alter the lumi- 
nosity functions, provided there are no populations 
which would be missing at any redshift. On the other 
hand, the level at which high-z populations are seen 
depends strongly on the assumed evolution of these 
galaxies. As an indication of the level of this effect, 
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we overplot predictions from the source count model 
of Rowan- Robinson (2001). The slight deviations of 
the model from the Euclidean slope is thought to be 
a numerical artifact (Rowan- Robinson priv. comm.). 



3.5 High-redshift evolution 

Can we use our local luminosity functions (section 

3.2) and the local spectral luminosity density (section 

3.3) to constrain the cosmic star formation history? 
An exhaustive analysis of the constraints is beyond 
the scope of this paper; however, we note here that the 
extragalactic background light is determined entirely 
by the evolving spectral luminosity density: 

where j v is the comoving spectral luminosity density, 
and Rodr the comoving distance element (Peacock 
2001; see also Hauser & Dwek 2001). A very large 
class of evolutionary models including pure luminos- 
ity, pure density, and mixed luminosity /density, sat- 
isfy 

j*{z)=j*(z = 0)f(z) (9) 

where f(z) is the Madau diagram normalised to 1 at 
2 = 0. We assume the far-infrared luminosity density 
is proportional to the volume-average star formation 
rate, following Rowan-Robinson 2001 (see equation 11 
below). A further assumption which underpins this, 
and indeed almost every study of the high- z sub-mm 
population, is that local galaxies can be found which 
are reliable templates of the high- 2 population. This 
is not necessarily the case, in which case equation 9 
will not hold. Nevertheless, if local galaxies are reli- 
able templates, equation 9 will hold in very general 
conditions. 

We can therefore use the observed extragalac- 
tic background light spectrum to constrain the dust- 
enshrouded cosmic star formation history, indepen- 
dently of constraints from redshift surveys of sub-mm 
blank-field surveys. This has been attempted by other 
authors (e.g. Gispert et al. 2000, Dwek et al. 1998) 
but without the benefit of a robust determination of 
the local luminosity density I v (z = 0). The authors 
therefore had to resort to assuming a small number 
of template SEDs, and explore the parameter space 
of f(z) allowed by range of SED models. We have 
the benefit of a very well-constrained determination 
of J„(z = 0). 

The determinations of the far-infrared extra- 
galactic background light historically have fluctuated 
in the literature, mainly reflecting the discovery and 
differing treatment of systematics (e.g. Puget et al. 
1996, Fixsen et al. 1998). We adopt the determina- 
tion of Lagache et al. (1999), since this is corrobo- 
rated by separate analysis of low-cirrus fields, and by 
the decomposition using the WHAM Ha survey and 
Leiden/Dwingeloo HI data (Lagache et al. 2000). We 



Figure 15. Example of the ±1<T constraints on the cos- 
mic star formation history defined in figure 13. Models are 
plotted which satisfy both the sub-mm source counts con- 
strain and the ±lcr likelihood range based on our analysis 
of the extragalactic background light. Note that the value 
of the peak redshift zt is strongly correlated with that of 
the evolution rate A, in that high values of zt require lower 
values of A, as can also clearly be seen in figure 13. Section 
4 discusses how these constraints be resolved with direct 
determinations. 



do not fit to DIRBE data shortward of 200/im, since 
these predictions rely mainly on an extrapolation of 
our SEDs in their rest-frames rather than an inter- 
polation. Warmer dust components are inevitable, re- 
gardless of starburst or AGN activity, but are not in- 
corporated into our SED models. We defer treatment 
of these warmer components, their associated colour- 
colour diagrams, and of the DIRBE points, to a later 
paper. 

We found that the extragalactic background light 
can be fit by a sequence of judiciously-placed 8- 
functions in f(z), and that there is no unique choice 
for these 5-functions. Therefore, a unique deconvolu- 
tion of f(z) from the extragalactic background light 
is not possible within the observational errors in the 
background. However, we can still make useful con- 
straints if we restrict the class of allowed forms for 
f(z). For simplicity, and to demonstrate the types of 
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Figure 14. Left figure shows the extragalactic background light (data from Lagache et al. 1999) modelled by the cosmic 
star formation history in equation 10. The data longward of 200/im plotted as broken lines is the FIRAS spectrum: full 
line is whole sky, and dashed line is Lockman Hole only. Also plotted are the DIRBE data points, not included in our 
fitting. The smooth curves are models of this spectrum, corresponding to cosmic star formation histories plotted in the 
right hand figure, provided our assumption about the evolving luminosity density (equation 9) holds. The full line is the 
global maximum likelihood fit quoted in table 3, and the long-dashed line has the same parameters except an enhanced 
high-2 star formation rate (B = 0.25) which is marginally (~ 2<r) inconsistent with the extragalactic background. The two 
short-dashed lines demonstrate selected alternative models: zt = 1 and zt = 2, both with A = 3 and B = 1. Note that in 
general the models with the higher star formation rates at z > 2 are also the models predicting the larger background at 
wavelengths A ~ 1mm (including in particular a comparison of the two short-dashed curves). In general, models with high 
volume-averaged star formation rates at z > 2 overpredict the sub-mm/mm-wave background. 



constraints which are available on /, we model the 
Madau diagram f(z) with the following simple pa- 
rameterisation: 



i.e., (1 + z) A power- law evolution to a threshold red- 
shift zt, and an arbitrary decline thereafter which we 
have chosen to be a straight line in the log f — z plane. 



f{z) = {l + z) 
f(z) = (l + z t ) A B z - 



! [z< z t ] 
[z > z t ] 



(10) 



The results of this fitting are shown in figure 
13, where we plot the confidence limits on A and z t 
marginalised over B. Restricting the parameter space 
explored to fixed values of B (i.e. instead of marginal- 
ising over all B values) favours different parts of the 
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Parameter Global peak 68% 



95% 



A 
B 



0.95 
4.20 
0.04 



0.75 1 



-+0.25 n 7c -+0.50 
-0.05 u -'°-0.15 

> 3.5 > 3.3 



< 0.52 



< 0.59 



Table 3. Maximum likelihood parameters for the extra- 
galactic background light fit. The parameters are defined 
in equation 10. The 68% and 95% confidence bounds 
quoted are in each case marginalised over the other two 
parameters with an 2 < A < 5 top-hat prior assumed. 
The global maximum is not affected by the prior on A. 



contoured region: low B values favour the low A / 
high z t region of the contours, while high B values 
favour the high A / low z t region. 

These constraints on the A and zt parameter 
space are determined only from the comparison of 
the extragalactic background light with our tightly- 
constrained determination of the local spectral lumi- 
nosity density, and made no use of (e.g.) the empir- 
ical constraints available on fi, or the 850/im source 
counts. These latter two constraints are overplotted in 
figure 13, for comparision with our extragalactic back- 
ground light constraint. The 850/im counts assume 
pure luminosity evolution to z = zt, in concordance 
with galaxy evolution at other wavelengths, and pure 
density evolution thereafter. O* is determined from 
the local 60pm luminosity density p6o with the fol- 
lowing relation (Rowan- Robinson 2001): 



P60 



to 



lL Mpc- 3 lOGYr 



(11) 



Some source count models of the sub-mm population 
have had difficulty in reproducing this fi„ constraint 
(e.g. Blain et al. 1999), though in our case we show 
below that there is good agreement with our extra- 
galactic background light modelling. 

Figure 13 can easily be translated into evolution 
tracks in the cosmic star formation history, as shown 
in figure 15. Note however that the parameters in fig- 
ure 13 are correlated, so that the error on the co- 
moving star formation rate at one redshift is corre- 
lated with that at another redshift, a fact which is 
rarely recognised (e.g. Gispert et al. 2000). Figure 13 
is therefore the "least-misleading" presentation of our 
constraints on the cosmic star formation history. 

An alternative to pure luminosity evolution 
sometimes used in the phenomenology of quasar evo- 
lution is to assume that the decline at the highest 
redshifts is pure density evolution, while still treating 
the increase from zero redshift as preserving pure lu- 
minosity evolution. We found this made only a modest 
difference to our 850/im 8mJy counts. These predic- 
tions are not plotted in figure 13 for clarity, but in 
the B = 0.04 and B — 0.6 cases, they are roughly 
equivalent to a la shift in the observed 8m Jy counts 
(i.e., a slight shift down and to the left of the 8m Jy 
counts constraint). The B — and B = 1 cases are 
non-evolving above the threshold redshift z t , so there 



are no changes in the 8mJy 850pm counts predictions 
in these cases. Also, the f2» predictions in general de- 
pend only on the evolution in luminosity density, and 
not on the comparative evolution of luminosities and 
number densities. 

The global peak of the likelihood in this parame- 
ter space is given in table 3, as are the 68% and 95% 
confidence bounds on each of the parameters when 
marginalised over the other two. We found accept- 
able fits with unphysically large values of A and very 
small values of z t (see figure 13) so we adopted an 
2 < A < 5 top-hat prior in calculating the confidence 
bounds. This did not affect the global peak. Figure 
14 shows the global best fit to the extragalactic back- 
ground light data. A notable requirement of these fits 
is a rapid decline in comoving star formation density 
above z t ~ 1. This is confirmed both by the fi* and 
850pm counts constraints, both of which are overpre- 
dicted in a B — 1 model (i.e. f(z) = constant above 
z = z t ). 

However, this star formation history is immedi- 
ately seen to be in conflict with determinations of the 
redshift distribution of sub-millimetre galaxies (e.g. 
Chapman et al. 2003) as well as other direct determi- 
nations. How can these be reconciled? We will address 
this in the following section. 



4 DISCUSSION AND CONCLUSIONS 

Our determinations of the local far-infrared to sub- 
mm luminosity functions span sufficient luminosity 
range to identify the convergence in the local spectral 
luminosity density. In particular, the 850pm local lu- 
minosity density is found to be 7.3 ± 0.2 x 10 19 hez, 
W Hz -1 Mpc~ 3 solving the "sub-mm Olbers' Para- 
dox" (Dunne et al. 2000). Our SED predictions rely 
on the extrapolation of the SLUGS colour-colour re- 
lations, both to higher and lower 60pm luminosities. 
Although we have no reason to suspect this extrap- 
olation is invalid, direct tests with sub-mm photom- 
etry of low and high 60pm galaxies are a key test 
of our methodology and are an attractive alternative 
to large photometric campaigns of local galaxies for 
constraining the local sub-mm luminosity functions. 

Using only the observed extragalactic back- 
ground light and our local spectral luminosity density, 
we derive a conditional constraint on the z < 2 ob- 
scured cosmic star formation history consistent with 
(1 + z) A evolution with A > 3.5 to z ~ 1. We find 
that the extragalactic background light, the 850pm 
8mJy source counts, and the fi* constraints all in- 
dependently point to a decline in the comoving star 
formation rate at z > 1. However, this is at variance 
with direct determinations, and with the constraints 
on the redshift distribution of sub-millimetre galaxies 
(e.g. Chapman et al. 2003). 

Perhaps the most uncertain of our assumptions 
is equation 9. This equation is true in a very gen- 



© 2002 RAS, MNRAS , 1-?? 



16 S. Serjeant & D. Harrison 



Figure 16. Predicted high-rcdshift cxtragalactic source counts at 170^tm and 850^tm for the A = 3.75, B = 0.04, zt = 1.0 
pure luminosity evolution model, compared to the observed counts from Barger et al. (1999), Dole et al. (2001), Scott et 
al. (2002), Cowie et al. (2002), and Webb et al. (2003). Note that none of these observed counts were used as constraints 
in generating the model, except the 850/im 8mjy count. In the 850/rni plot, each of the data points are from a different 
source and (unusually for a cumulative source count plot) are therefore statistically independent, with the exception of the 
> 8m Jy data which is all from Scott et al. (2002). The overprediction of the faint 850/im counts is further indication that 
the evolution at faint luminosities is less strong that that at bright luminosities (section 4). 



eral class of cases, including pure luminosity evolu- 
tion, pure density evolution and any mixed luminos- 
ity/density evolution models (i.e. those in which the 
shape of the luminosity function is invariant with red- 
shift, but not the break luminosity or the number 
density at the break). However, if (for example) the 
high-luminosity sub-millimetre galaxies evolve more 
quickly than their low-luminosity counterparts, then 
equation 9 will not hold. Stronger evolution at the 
high-luminosity end would give the comoving lumi- 
nosity density a warmer colour temperature at high 
redshift, which in turn would permit higher numbers 
of sub- millimetre galaxies at high redshift without vi- 



olating the extragalactic background light constraint. 
The warmer colour temperature may also explain why 
our models fall short of the data at 140/im (figure 
14), though this may also be due to our redshifted 
SEDs being less well-constrained at the shortest wave- 
lengths. The f2» constraint would also be alleviated, 
since the extrapolation to galaxies with lower lumi- 
nosities would yield less additional star formation. 

A further indication that the evolution may 
be stronger at high luminosities is provided by the 
predicted blank-field high-redshift source counts at 
170/im and 850/im, shown in figure 16. It is im- 
portant to stress that none of the observed source 
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counts in this figure were used as a constraint on the 
model, except the 850/im 8m Jy count. While the ob- 
served 170/im and bright 850/im counts are both en- 
couragingly well-reproduced by the model, the faint 
850/im source counts are over-predicted. Suppressing 
the strength of evolution in the less luminous popu- 
lation would have the effect of reducing the predicted 
850/im counts at the faint end. 

An alternative solution to alleviating the con- 
straints on f2* is that the initial mass function is 
more top-heavy at high redshifts. The infrared lu- 
minosities are dominated by > 6Mq stars (Condon 
1992), so altering the initial mass function can have 
a strong effect on the fl, calculation while leaving 
the extragalactic background and source count pre- 
dictions roughly intact. There are good prima facie 
reasons for supposing this might be the case (e.g. Lar- 
son 1998 and references therein), though if it were to 
entirely explain the discrepancy with the 850pm red- 
shift distribution, the top heavy initial mass function 
would have to exist in low luminosity high-redshift 
galaxies, as well as ultraluminous high-redshift galax- 
ies. 

Within the cosmic star formation histories f(z) 
we have considered, we can therefore conclude that 
the initial mass function is very different at high red- 
shift, and/or some differential evolution is required; 
all pure luminosity, pure density and mixed lumi- 
nosity/density evolution models are excluded by the 
850/im blank-field survey redshift distributions if the 
initial mass function is the same at high redshift. Fur- 
thermore, comparison with the faint 850/im source 
counts already indicates the presence of differential 
evolution. 

The issue of mm-wave/sub-mm/radio photomet- 
ric redshifts is very timely given the advent of the 
SHADES survey (van Kampen et al. in prep., Mortier, 
Serjeant et al. in prep.), which will combine SCUBA 
and BLAST data with the aim of obtaining photo- 
metric redshifts accurate to Az ~ 0.5. The accuracies 
which are obtainable are still the subject of debate 
in the literature (e.g. Hughes et al. 2002, Blain et 
al. 2003) and are limited by the availability of local 
templates. It would be interesting to address the pho- 
tometric redshift accuracy from the point of view of 
our SED ensemble. However, our results have already 
indicated that local ultraluminous galaxies may not 
necessarily be templates for the high-z ultraluminous 
systems being found by sub-mm blank-field surveys. 
The photometric redshifts will therefore need careful 
calibration against spectroscopic redshifts (e.g. Chap- 
man et al. 2003), obtained for example using Spitzer 
identifications (e.g. Egami et al. 2004, Serjeant et al. 
2004, Ivison et al. 2004, Frayer et al. 2004, Charman- 
daris et al. 2004). 
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